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We study, analytically and numerically, phase locking of driven vortex lattices in fully-frustrated 
^ . Josephson junction arrays at zero temperature. We consider the case when an ac current is applied 

O ' perpendicular to a dc current. We observe phase locking, steps in the current- voltage characteristics, 

with a dependence on external ac-drive amplitude and frequency qualitatively different from the 
5^ , Shapiro steps, observed when the ac and dc currents are applied in parallel. Further, the critical 

CIh' current increases with increasing transverse ac-drive amplitude, while it decreases for longitudinal 

^ ' ac-drive. The critical current and the phase-locked current step width, increase quadratically with 

^ , (small) amplitudes of the ac-drive. For larger amplitudes of the transverse ac-signal, we find windows 

where the critical current is hysteretic, and windows where phase locking is suppressed due to 
dynamical instabilities. We characterize the dynamical states around the phase-locking interference 
condition in the IV curve with voltage noise, Lyapunov exponents and Poincare sections. We find 
that zero temperature phase-locking behavior in large fully frustrated arrays is well described by an 
effective four plaquette model. 



PACS numbers: PACS numbers: 74.81.Fa, 74.25.Sv, 74.25. Qt 



I. INTRODUCTION 

> 

' Phase locking phenomena are found in a wide variety of nonlinear driven systems in condensed matter physiesi 
It takes place when an internal frequency of the system locks to a rational multiple of the frequency of an external 
ac-drive. A simple example of this is the case of an overdamped particle moving in a tilted washboard potential, 
where the frequency of motion of the particle over the periodic potential can be locked to multiples of the frequency 
of a superimposed ac force for a finite range of the dc force (tilt of the washboard) . Since the internal (washboard) 
frequency is proportional to the mean velocity of the particle, phase locking results in a constant mean velocity for a 
certain range of dc-force curve when the interference condition is satisfied. A particularly well known realization of this 
effect is Shapiro steps^ in the dc current- voltage {IV) characteristics of a single small area Josephson junction driven 
by a time periodic current. Within the washboard analogy outlined above, a simple analysis provides expressions for 
'"^ ' the appearance of Shapiro steps at specific voltages^ corresponding to integer multiples of the driving frequency. 

' Driven systems with many degrees of freedom can also exhibit phase locking. This has attracted broad scientific and 
O technological interest since phase-locking in complex systems can either be induced by collective effects, providing 
for a low dimensional interpretation of the phenomenon, or itself induce collective (low dimensional) behavior in 
the complex system. Phase-locking experiments have provided information about such dynamical response of non- 
equilibrium collective states, where dimensionality, thermal fluctuations, quenched disorder, and the magnitude of 
external fields can be very relevant. A particular well known example is the large Josephson junction array (JJA), 
' with N X N junctions, driven by an external current (/^c + IacCOsflt)x with frequency Q and with an applied 
magnetic field density / = Ha^/^o, where H is the magnetic field, a the lattice period of the Josephson array and 
$0 the quantum of flux. Giant Shapiro steps at voltages — NnKl/2e, n being an integer, have been observed 
experimentally in zero magnetic field (/ = 0)^ Fractional giant Shapiro steps at voltages Vn,q = Nnhn,/2eq, were 
observed experimentallj^*S for strongly commensurate magnetic fields, / = p/q, with p, q integers, and extensively 
investigated in numerical simulationsiSsSiifliii Also, subharmonic giant Shapiro steps at voltages Vn^m — Nnh£l/2em 
were observed experimentallyS. for zero magnetic field (/ = 0) , and attributed to the nucleation of complex collective 
dynamical statesiSiH induced by disorder or inductance effects. Shapiro-like phase-locking is also observed in the case 
of driven vortex lattices in bulk superconductors with two-dimensional periodic pinning arrays, as recently reported 
both experimentally-'^'* and theoretically.*^ Also superconductors, where vortices are driven over a one-dimensional 
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potential generated by thickness modulationsji* or are confined to move through mesoscopic channelsjii show Shapiro- 
hke phase-locking. Moreover, systems with many degrees of freedom in the presence of quenched disorder can also 
exhibit phase-locking when there is a dynamically induced periodicity, like charge density wavesi^ and vortex lattices 
in superconductors with random pinningii2iSS*2i 

In the phase-locking examples mentioned above, the ac-drive is applied parallel to the dc-drive. However, it was 
recently shown that a different type of phase-locking, distinct from Shapiro phase-locking, is possible in vortex lattices 
if the ac-force is applied perpendicular to the dc-force.^^-^'^ In this case the interference effect is due to an effective 
parametric ac drive in the longitudinal direction, which is induced by the transverse ac drive. In several systems, 
like charge density waves or single degree of freedom systems {e.g., the single Josephson junction), the dynamical 
variables are such that perturbations or displacements can be induced in only one direction {i.e., the displacement 
field is a scalar). An important characteristic of vortex lattices in superconductors is that the displacement field is two 
dimensional. In particular, the behavior of displacements in the direction perpendicular to the driving force shows 
phenomena like a transverse critical currenlj^iSiSSiSi and a transverse freezing transitio n^^'^^ at high velocities. Phase 
locking in ac-driven vortex lattices, where the ac and dc forces are perpendicular, arises as a direct consequence of 
the nonlinear coupling between the two directions of motions. 

Transverse phase-locking has been reported for vortex lattices moving in rectangular or triangular pinning arraysSS 
as well as in arrays of randomly distributed pinning centers. ^'^ In this paper we investigate the possibility of transverse 
phase locking in a two dimensional (2D) fully frustrated JJA, where the average of the external magnetic field 
corresponds to one half flux quantum per plaquette, / = 1/2. This system has several attractive properties. The 
presence of a magnetic field (/ 0) breaks the axial symmetry in the direction of the bias current, and 2D-cooperative 
behavior may come into play. This leads to the well-known fractional giant Shapiro gteps^iSiL^iSiiSiii induced by a 
longitudinal ac current. It also, as we will demonstrate in this paper, allows for transverse phase locking when the 
ac current is perpendicular to the dc current, since the two directions of motion become coupled. Non-equilibrium 
dynamical phases for fully frustrated JJAs driven by a dc current have previously been studied»22i^ Phase locking 
can be used to characterize tem poral order in the different dynamical phases of the JJA at high velocities by their 
ac-response, as was done in Ref. 20,21 for bulk superconductors. 

Here, we report transverse phase locking steps in the IV characteristics, similar to the longitudinal giant Shapiro 
steps, but with very different characteristic dependencies on external ac-drive amplitude lac and frequency fi. The 
critical depinning current in the system with transverse ac drive is larger than the critical current of the dc driven 
system. For lac/^ ^ 1 the depinning critical current Ic and the phase-locked step width A^i for V = ilh/2e increase 
quadratically with lac- For lac/^ > 1 we find windows of lac/^ where depinning is hysteretic and the periodic phase- 
locked dynamical states become unstable. We characterize the dynamical states around the phase-locking interference 
condition in the IV curve with the voltage noise, Lyapunov exponents and Poincare sections. We find that zero 
temperature phase-locking behavior in large fully frustrated arrays is well described by an effective four plaquette 
model. 

The remainder of this paper is outlined as follows. In section II we present the model used for simulating the 
dynamics of the fully frustrated JJA. In section III we develop an analytical framework for predicting critical current 
and phase-locking properties of the fully frustrated JJA. Section IV presents simulated transverse phase-locking steps 
in typical IV curves obtained from an effective four plaquette model for the JJA. We calculate numerically the 
dependence of the critical current and the magnitude of phase-locking steps with the amplitude lac and frequency 
n of the external ac-drive. The results obtained are analyzed in more detail by studying voltage noise, Lyapunov 
exponents and Poincare sections for the dynamical states around the phase-locking interference condition. We also 
compare numerical simulations results for large arrays with those obtained using the effective four plaquette model. 
The discussions and conclusions of the investigation are presented in sections IV and V. 



II. MODEL 



We study a current driven JJA with an ac current perpendicular to the dc current, as shown in figure^. A magnetic 
field, H, is applied such that half a flux quantum, / = Ha^/^o = 1/2 with being the area of a plaquette and 
$0 — h/2e being the flux quantum, penetrates each plaquette; corresponding to the fully frustrated XY modelr^^ 
where the ground state is a "checkerboard" vortex lattice, in which a vortex (flux quantum) penetrates every other 
square grid (see Fig. ^). In such ground state, current and phase differences in the junctions are described by a 
repeated two-junction by two-junction (2x2 plaquette) superlattice unit cell. 

Numerical simulations^'* of large driven arrays suggest that this spatial periodicity is preserved when the dynamics 
is phase-locked to an external ac-perturbation applied in parallel to the dc force. We will show later, in section IV, 
that this is also a good approximation when the ac current is perpendicular to the dc current. We will therefore 
consider the simple system of a 2 x 2 superlattice unit cell of the array and the associated gauge-invariant phase 
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differences in a field of / = 1/2, with the dc current (per plaquette) Idc perpendicular to an ac current (per plaquette) 
/acC0s(r2t), as shown in Fig. |21 Flux quantization, total current conservation at the central node, and the applied 
total currents in the two directions give the following governing equations, 

a + /3-(5-7 = 7r(l + 2?i) (1) 

/3 + 7 + sin/3 + sin7 2IacC0snt (2) 

d + (5 + sin a + sin 5 = 2Idc (3) 

d + 7 — /3 — (5 + sina + sin 7 — sin /? — sin (5 = 0, (4) 

where n is an integer, t is the normalized time in units of to — 2eIoRN/h, i?jv being the normal state single junction 
resistance, 51 is the normalized frequency in units of l/to, lac and Idc are the normalized external currents in units 
of the single junction critical current Iq. This model was introduced by Benz et al.'^^ (for lac — 0) to study the dc 
current- voltage curve of the f — 1/2 array. They obtained analytically that the critical current per junction of this 
model is 4 = {y/2-l)IoM The same model was later used in Ref. 'o'to study the (longitudinal) Shapiro steps. Since 
the analysis done in the work of Refs. i9ii33i did not include the transverse ac current with accompanying transverse 
voltage drop, an additional constraint of P = 6 was implied, reducing the model system to two dynamical degrees of 
freedom. In contrast, our model system of a four plaquette unit cell consists of three effective dynamical variables. 
We calculate the instantaneous longitudinal and transverse Vy (normalized) voltages per junction as, 

V, = ia + 5)/2 (5) 
Vy = (/3 + 7)/2, (6) 

and the IV characteristics, Vx = (Vx) as a function of Ida where (••■) is a time average. The total longitudinal 
voltage vt for an. N x N array, built with this 2x2 superlattice unit cell, is vt = Nvx- When vortices move with 
a mean velocity u in such 2a x 2a superlattice structure, we can obtain the normalized voltage Vx using the relation 
2'Ku/2a = Vx, where a is the array periodicity (see Fig. ^). The intrinsic washboard frequency for vortices moving 
with velocity u in the periodic potential of the JJA is luq ~ 2TTu/a. Phase-locking in the longitudinal direction is 
obtained when the frequency of the ac drive locks to a rational multiple of the intrinsic washboard frequency. For the 
71-th harmonic this corresponds to luq — nQ, i.e. 2TTu/a — nQ. This leads to phase-locking at voltages Vn,2 = {n/2)Q 
for fully frustrated JJA. In general, for / — p/q, the ground state has qa x qa superlattice structure, therefore the 
voltage for vortices moving with velocity u is 2'Ku/qa — Vx, and phase- locking for the n-th harmonic is obtained at 
voltages Vn^q = {n/q)n. This is the condition for the so-called "fractional giant Shapiro steps"i^*S*2i2i2iiSiii. 

III. PHASE-LOCKING AND CRITICAL CURRENT ANALYSIS 

We will in this section assume that the dynamics of the system is represented by the simple two-plaquette degrees of 
freedom as shown in figure|21 We will apply the following linear transformation of the phase- variables of Eqs. C3l-I01l: 

= ^ (7) 

= ^ (8) 
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= (9) 

= ^ • (10) 

With these variables we can represent the constraint [J — ^) oi Eq. ^ as x + "^y ~ and thereby write the 
relevant three degrees of freedom in either of the two following forms, eliminating '^x'- 

<i>:r + sin f y sin $^ Idc (11) 

^y + cos sin $,y = I ac COS fit (12) 

2i'y - cos^xcos'^y + cos^ySin^y = 0, (13) 
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or eliminating 



cos r sin 



+ sin 5*2^ sin = cos fii 



2i'x + cos sin " cos cos = 



(14) 
(15) 
(16) 



The normalized voltages are given by, Vx — and Vy = ^y. We will in the spirit of the usual Shapiro analysis 
assume that 



— —— sin Qt 

which is the solution to Eqs. H12|l and p5|l for large lac and f2. 



(17) 



Critical Current 



We will here look at Eqs. ifT^ and lfTB|) . Assuming <i>x = $1"^ and 
constants and |e(t)| ^ 1, the static contribution to equation ((TE|l is: 



cos$i°'sinvI/W = Jo(^:^ 
^ tanvI/W = 



01 -^^cos^-j,") 



cos$i°^ 



^"^ '-e(i), where $1"^ and ^-l,"^ are 



(18) 
(19) 



where J„ is the n'th order Bessel function of the first kind. Inserting this into the static part of equation (|14|l yields 



Idc 



cos < tan 



cos $ 



(0) 



sm 



$(0) 



(20) 



This expression provides a unique relationship between the constant phase, <I>i°\ and the dc current, Idc- However, 
for increasing Idc, there exists a critical value, Jj, for which no real <I>i°^ can satisfy Eq. J^O)). This value is given by 



/I = 



cos < tan 



cos$ 



(0) 



sm 



$(0) 



(21) 



which is the predicted critical dc current for static states {vx = {^x) = 0). We notice that the identical expression for 
the critical current, Eq. H21() . can be obtained from the anisotropic dc driven system. 



a + /3 - (5 - 7 = 7r(l + 2n) 

/3 + 7 + rsin/3 + rsin7 = 

a + 6 + sin a + sin 6 = 2Idc 

d + 7 — /3 — (5 + sin a + r sin 7 — T sin /3 — sin 6 — 



(22) 
(23) 
(24) 
(25) 



where the anisotropy, T (suppression of transverse critical current), is given by the standard Shapiro^ critical current, 



r = ^o(%) 



B. Phase-Loclting 

We will here use equations (|ll() - (|13|l . since {'i'y) = 0(mod7r) for {<^x) 7^ provides for a simple description of the 



dynamics. We will assume the ansatz, $ 
reads: 



(0) 



2* 



rtt and 'i'y = A sin fit + B cos fit. The equation for 'i'y, H13|) . now 
t/y ~ cos($(°) + ^t) cos ^y + sin 'J'y ^ Jfc ( j cos fcf7i = , (26) 
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where we will use the approximations: cos^'j, « Jq (V^^ + i?^) and sinVt^ « ^y. With the ansatz for above, this 
equation has no static component. The time varying component, at frequency f2, yields the coefRcients A and B 

21] cos $1*'^ - [Jo (iff) + J2 (iff)] sin$l"^ 



A, = Jo ( ^Al + Bi] ^ ;VK.:^ [ V^o + Bi ) (27) 



+ Jl (if) - j| (^1 



. .oJo(v^^.-o^) ^ ^"^^-t:,^1;%^r1iM^^^^" Mv^^^^^) ^ (-) 



where and Bq are the solutions for Jq (^■\/Aq + _Bq j = 1. Thus, the solution [A,B) — (AojSo) is correct up to 

order fi"^ and {A,B) = {Ai,Bi) is correct up to . Inserting this solution {Ai,Bi) for into the ^'^-linearized 
equation l(TT|) gives the static properties: 

Jn ( JaI + BI^ i (Ao cos + Bo sin = /dc - ^ (29) 

— r\ 1 217- J2 (%) sin2$i°^ 1 

= Jdc-f^. (30) 

Thus, the locking range for the this step can be found to second order in l^j^l- The dominant part of this expression 
for the range in phase-locking yields: 



A5i 



^^('^) 31]^ + l(Jg(^)+J|(^)A _ (3^^ 



+ 4 (^) - j| (^) V {m^ + 4 (%) - 4 (%)) 



The expression displays quadratic growth of the phase- locked step size for small lac- This is consistent with the particle 
(pancake) model results^^ for vortices in rectangular pinning arrays, and it is different from the known longitudinal 
(Shapiro) phase-locking of JJAsA^iSil 

In addition to the range of phase-locking, AS'i, equation (|30|l provides information about the offset, A/i, of the 
phase- locked step relative to the Ohmic (linear) curve (see inset in figure OJ. The offset is given by the part of the 
equation that does not depend on '^^x^. From Eq. (|3(J|) we have: 

Ar = » (, „ 1 417^ + Jg(^)+2J|(^) \ 

Expressions (31) and (32) provide a second order (in 12"^) description of the phase-locking step magnitude and 
location as a function of the system parameters, 17 and lac-, for large 17. 



IV. RESULTS AND DISCUSSION 



We will here show the results of numerical simulations of the system analyzed in the previous section. The sim- 
ulations are conducted with numerical parameters corresponding to the model parameters, using a fourth order 
Runge-Kutta method such that the normalized time step typically is no larger than 1% of the period of the driving 
frequency, and often smaller. Since we are mostly concerned with dc lY characteristics, we choose to acquire data 
for averaging over many ac-periods of motion (typically 10^-10'^) after a sufficient initial time of interval allowed for 
transient behavior. Simulated IV characteristics are obtained by performing the necessary averages as described, and 
then changing the dc current, Idc-, slightly to acquire the next point on the IV curve. All simulations are conducted 
for the fully frustrated case of / = \- 

Figure O shows a simulated IV characteristic, v^, as a function of /^c, simulated for 17 = 1 and lac — 2.3, for the 
simple four plaquette model showed in figure 12 described by the equations is obvious from the figure, 

we obtain clear signatures of critical current(s) and phase-locked steps. Specifically, we observe the AS'i step at 
Vx = 17, and steps AS'i and AS'2 at Vx ~ ^17 and Vx — 217, respectively. We observe a critical current larger than the 

previously predictec^ value of Ic — — 1 w 0.41 for a fully frustrated dc-driven system. The characteristics of this 
plot are very similar to the behavior observed in JJA with parallel ac -I- dc drives, obtained both by simulations^iifl 
and experiments^, as well as analytically for the four plaquette modeli^ 



6 



Comparing simulations of the fuU-RSJ dynamical equationsSLSS for different large arrays {N x N junctions) we will 
later (below) demonstrate that the simple model of a 2 x 2 array gives very good description of the JJA dynamics. 
However, we will first compare the predictions of the analytical treatment of the previous section to numerical results. 



A. Critical Current and Phase-Locking for > 1 



In order to verify the simple theory for critical current and phase-locking behavior developed above, we have 
conducted numerical simulations of the four plaquette system described by Eqs. Q-Q for 1 < 17. 

The first set of simulations are conducted to investigate the critical current, Ic, as it is described by Eq. H21|l . 
This expression provides an estimate of the critical current, /J, for which the JJA switches from a zero- voltage state 
((14) = and Idc < Ic) to a non-zero voltage state. The simulations are conducted accordingly, starting the system at 
rest for small Idc and slowly increasing the dc-bias until non-zero average voltage is detected. The results are shown in 
figure 0^, where the solid curve represents the expression, Eq. H21|l . and the markers represent the simulation results 
for several frequencies 1 < < 3 as a function of the characteristic ratio. Iff. The size of the markers are larger than 
the error on the estimated critical current. It is obvious that the agreement is very good for all simulated data sets, 
and we conclude that the critical current, as given by Eq. (|21|l . is a relevant estimate for not smaller than 1. 

Figure^ shows the critical current, evaluated from numerical simulations when the dynamical system switches 
from the non-zero voltage state to the zero- voltage state (see inset in figure O)). We have here shown the results of 
numerical simulations with markers as for figure^^, together with the solid curve of Eq. H21|l . However, it is clear from 
the figure that the critical current, for decreasing dc-bias may be smaller than for increasing dc-bias (/J > /J). 
Since this is a multi-dimensional system, the critical current may be hysteretic, such that decreasing the dc-current, 
Idc, for non-zero voltage states {v^ — ($3,) ^ 0) is subject to different critical characteristics. One simple way of 
investigating this is to assume a non-phase-locked state of voltage Vx 0, such that = Vxt- A primitive analysis 
can provide a hint to this hysteresis. 

The critical current analysis of the previous section, is obviously a critical current for a system operated at the 
(x) = branch of the IV-cuive. We may instead analyze what may happen for a non-phase-locked {^x) = ^'x 7^ 
state. We will still assume (17) to be an appropriate description of the transverse current. However, equation H13f) 
becomes (for small and with no resonance to fl): 



2*. 



COSVxt -\- Jq 



j2(^)+4z;- 



■ cos Vxt 



72 (lao 



■ siniiri . 



(33) 
(34) 



Inserting the ^'y solution into equation yields the static component: 



Jl 



2Vrr 



(35) 



where the second term on the left hand side is the result of the resonant mixing between the propagation, {^x) = Vx 
and the transverse oscillation, ^f^. However, the overdamped dc driven pendulum equation is also subject to the 
following simple relationship 



+ (/i)2 - 

Combining the two expressions provides the relationship 



Idc ■ 



(36) 



Jl 



2v,. 



72 

■^0 



< 0.826591 . 



(37) 



As we have indicated, the critical current, /J- has a maximum value at around 0.82 (for Jq (%^) = and ~ 0.33). 
Thus, we can argue that propagating {{^x) 7^ 0) solutions may exist for Idc > Ic ~ 0.82, which provides for a 
hysteretic IV characteristic switching between zero and non-zero voltages in the range I^ < Idc < Ic 1 when I^ < ■ 
Notice that when Jj- > /J , the relevant critical current for both zero and non-zero voltage states must be /J , since no 
static states exist for Idc > Ic- However, when I^ < l}, the actual critical current for switching into a zero- voltage 
state may be anywhere in the interval [/^; /J]. 
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Figure 0] clearly indicates the hysteretic switching in the IV characteristics when < Idc < ll , which is the case 
for smaU |Jo(%^)|- 

We note that the above rather primitive analysis of the hysteresis provides a fairly good agreement with the results 
of numerical simulations. The results of the analysis are not completely consistent with its assumptions in that the 
resulting amplitude of for the optimized Vx ~ 0.33 is about 1.5, which is not a small number. However, a more 
detailed (nonlinear) analysis of a single frequency representation of '^y yields quantitatively similar and qualitatively 
identical results (/J < 0.77) as the above, and we therefore conclude that the simple explanation for hysteresis 
presented here is relevant. 

For Jo (%^) = we can provide an explicit approximate expression for by assuming the critical current is given 

by the value of Vx which optimizes Ji (^^it) • This leads to Vx ~ 1/3.6, which inserted into the above equation yields 

the optimized coordinates: {I^,Vx) = (0.81,0.28). 

It is noteworthy that we observe, as predicted by the expression for /J, Ic{i^,Iac) to be larger tha.n the dc-driven 
system, Id^^ lac) > -^c(^^, 0) ~ 0.41. Hence, a transverse ac driving leads to an enhancement of the critical current. 
This is contrary to the case with the ac-current parallel to the dc-current, where the critical current is reduced, i.e. 
Ic < lc{^,0) « O.4I.W1I" 

The predicted range, A^i, in Idc of phase- locking, as given by Eq. 1)31(1 . is investigated through simulations similar 
to the above study of the critical current. Comparisons between the predicted expression and results of numerical 
simulations are shown in figure [5^, which displays the largest magnitude of the range in dc current for which phase- 
locking is observed as a function of the characteristic ratio, for different values of Q in the interval 1 < < 3. 
Markers represent results of numerical simulations and solid curves represent the corresponding predicted results of 
Eq. H31|) . It is obvious that the simulated parameter sets provide very good overall validation of the perturbation 
analysis, with the larger of the simulated frequencies providing better agreement than the smaller, as expected. 
However, an observation common to all simulated frequencies is that large deviations from the expected behavior are 
found for parameter values {Q and lac) leading to Jo{^) ~ 0. The reason for this discrepancy is likely due to a 
dynamical instability, which can be explained by the perturbation analysis above. The average equilibrium position, 

{^y), of the variable '^y can be observed from Eq. H26|l if we write '^y — ^y"'' -I- ipy, where ^'y"'' is varying slowly in 
time (much slower than fi), and ipy represents all high frequency (including fl) contributions (|V'j/| ^ !)• The slow 
evolution of equation \2<i\i can then be written, 

2^i,°)+Jo(^)sinvI/(") = 0. (38) 

Thus, we find that the stable position of 'i'y is: 

The consequence of this abrupt transition in ^'^ is that the locking phase, must experience a similar abrupt 

transition of tt, as can be seen from equation ((ll|l . We therefore claim that the apparent discrepancy observed between 
the numerical simulations and the perturbation theory near the roots of Jq is a result of dynamical instabilities 

arising from switching the average phase, {^y) between and tt. 

We finally show the comparisons of the center of the phase-locked step as a function of the characteristic ratio, 
for different values of fl in the interval 1 < Q < 3. The predicted behavior, Eq. (|32|l . is subject to the same 
issues as the predicted range of phase- locking since the origin of both expressions is Eq. H30|) . Figure \^ shows the 
offset, A/i = /i — ri, between the center of the step, Ji, and the Ohmic curve. As is the case for the phase- locking 
range shown in figure the comparison between numerical simulations (markers) and the corresponding predicted 
offsets (solid curves) is very good, except for the instabilities near the roots of Jq (%^)- Notice that the phase- locking 
analysis leading to the predictions Eqs. (|31|1 and (|32|l does not depend on the sharp transition between and 
TT. The reason is that this transition provides only a sign change in the effective equations of phase-locking, and the 
magnitudes of locking range and offset are therefore unaffected as long as Jo (%^) 7^ 0. 

Based on the above presented comparisons between the numerical simulations of critical currents, range of phase- 
locking and position of the phase-locked step in the IV characteristics of a transversely ac-driven JJA and the 
corresponding results from simple perturbation analysis, we conclude that the high frequency behavior is well described 
by the presented analytical treatment. 
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B. Critical Current and Phase-Locking for intermediate and low 

In Fig. 121 we show the critical current behavior for intermediate and low frequencies. For intermediate frequencies 
(figure IHK) we observe how the critical current, /J, increasingly deviates from the high frequency behavior outlined 
above. Even so, we notice that the overall behavior of the critical current is qualitatively well described by the analysis 
leading to equation (|21() for > i. We have, for comparison, included an example of the critical current for the 
longitudinally ac-driven JJA as an inset. Not surprisingly, decreasing the frequency further (see figure EId) results in 
rather large discrepancy between the high frequency analysis of section IIIA and the numerical simulations, and no 
universal behavior of the critical current as a function of the characteristic ratio, can be found. However, we do 
observe that the critical current does seem to increase quadratically for small lac- 

In Fig. Owe show the phase-locking range, AS*!, at = as a function of lac/^ for intermediate (fi > 0.5) and 
low frequencies (f2 < 0.5). Again, as for the critical current we observe that the intermediate frequency range provides 
for reasonably good qualitative comparisons between numerical simulations and the high frequency analysis of section 
IIIA. We have, for comparison, included an example of the comparable range of phase-locking for the longitudinally 
ac-driven JJA as an inset. A noticeable feature of figure Et- is that the dynamical instability discussed above around 
Jo (■^) =0 seems to widen as the frequency is lowered. Figure [TJ) shows how this instability provides for increasing 
discrepancy between high frequency analysis and numerical simulations. However, we notice that even the very low 
frequencies retain the basic feature of quadratic growth of the phase-locking range as a function of lac for small lac- 

C. Dynamics of Phase-Locking 

Let us now analyze in detail the dynamics of the voltage responses (increasing and decreasing dc-current) to elaborate 
on our previous results: critical current hysteresis and windows without transverse phase locking. We calculate IV 
curves and Lyapunov exponents as a function of lac and f2. In order to distinguish between periodic or quasi-periodic 
dynamics and chaotic dynamics we calculate the maximum Lyapunov exponent. A, following the standard methods 
of nonlinear dynamics , This means that a small perturbation, e(0), to the initial condition will displace the new 
trajectory by an amount \e{t)\ ~ |e(0)|e'*'*. The Lyapunov exponent is then defined as 

A = lim i In = hm \(t) . 

To recognize a chaotic trajectory we evaluate the maximum Lyapunov exponent. If A > the trajectory is locally 
unstable; i.e., initial points that are arbitrarily close to each other are macroscopically separated by the flow after a 
sufficiently long time and the attractor is chaotic. Negative Lyapunov exponents are obtained when trajectories that 
start sufficiently close to a subset are attracted to it. Here we will show two particular cases: lac/^ = 3.0/1.5 = 2, 
corresponding to a set of parameters where no hysteresis is obtained, and lac/^ = 4.05/1.5 = 2.7, corresponding to 
the hysteretic regime. In Fig. |Hlwe plot the IV curves and maximum Lyapunov exponents for fl = 1.5 and lac = 3.0 
(a-b) and lac = 4.05 (c-d). The exponents are estimated from A w X{t) after a finite time t = 1024r, with T = 2tt /Q,. 
For lac = 3.0 we show a range in I^c where a wide transverse phase locking step exists (Fig.|S^), and the corresponding 
maximum Lyapunov exponent. A, is shown in Fig. (SJa. We see that within the step we have A < 0, with the most 
negative value of A at the center of the step. Outside the steps, the Lyapunov exponent is nearly equal to zero, 
A < 0, corresponding to quasiperiodic behavior. A different behavior is obtained for lac = 4.05, shown in Fig. (S];, 
where we see that the step disappears and thus there is no transverse phase locking in the same Idc range where we 
find a step in Fig. The maximum Lyapunov exponent, plotted in Fig. is small but positive for the Idc range 
around the expected location of the step, the smallness of A implies that the dynamics can be either chaotic (A > 0) 
or quasiperiodic (A = 0) in this case. The IV curves near the corresponding critical currents are shown as insets. 
We see that the absence of hysteresis in critical current is associated with the occurrence of transverse phase locking. 
Inversely, hysteresis in critical current is obtained for approximately the same parameters for which transverse phase 
locking is absent. This is in agreement with the above analysis that indicates the critical current hysteresis is present 
in the vicinity of Jq (■^) =0, which is also the location of the dynamical instabilities of the locked phase of the A5'i 
step. 

In summary, around the transverse phase locking step region we can distinguish three different voltage responses: 
A, B and C, which are indicated in Fig. |S1 We now calculate the voltage power spectrum and Poincare sections 
to distinguish these three types of dynamical behavior. This way to characterize dynamical behaviors was used 
before in capacitive rf-biased JJAi^SiS We analyze both transverse and longitudinal voltage power spectra. From the 
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instantaneous transverse voltage we obtain the transverse voltage power spectrum: 



1 f'^' 2 
'S'y(^) = 7f7 / dtVy{t)exp{iujt) 



(40) 



t Jo 



where Tt — NtA.t. From the instantaneous longitudinal voltage, Vx, we obtain the longitudinal voltage power spectrum: 



For studying the nature of the attractor in the different regimes it is useful to consider a Poincare section of the 
phase-space trajectoriesi^ We consider the stroboscopic Poincare section of the trajectories in the d^x/dt vs. sin<i>a; 
plane, recording the values taken by these variables each period of the ac drive. In Fig. I^we show the power spectra 
and Poincare sections for = 1.5 and for the lac and I^c values corresponding to the A, B and C regimes. For 
each case we show the longitudinal, Sx{(^), and transverse, Sy{uj), voltage power spectra as a function of Lu/ft and 
their corresponding Poincare sections. Let us first discuss the case corresponding to the regime B, in which there is 
transverse phase-locking. This is shown in Fig. ^jp for Idc = 1-66 and lac = 3.0, which corresponds to the step with 
mean voltage — (Vx) — (see Fig.|SJi, regime B). We see that the longitudinal power spectrum, Sx{u!), presents a 
delta-like peak for uj = 2n. Thereby, the first harmonic of longitudinal voltage fluctuations is locked to 20., as expected 
for this step, since it corresponds to n = 2 in the phase-locking condition loq — nfi. The phase-locking with a double 
frequency corresponds to the case when the vortex lattice oscillates in full synchrony with the transverse ac-current, 
and the ground state repeats itself after one period of the ac drive. In the transverse voltage power spectrum, Sy, 
there is a sharp peak at 11. This is characteristic of transverse phase-locking: the dynamics in the transverse direction 
locks at half the frequency than the dynamics of the longitudinal direction. This is so because in a single period of the 
ac-drive, T — 27: /Q., the longitudinal component moves forward n steps in the lattice period a, while the transverse 
component completes only the first half of its oscillation. In Fig. we show the Poincare section corresponding to 
this case in the regime B. The figure shows a very localized Poincare section since the trajectory always comes back 
approximately to the same location in phase space in each ac cycle, since the trajectory is periodic (closed orbit). 

Now we analyze the case corresponding to the A regime, which is for a current outside the step, Idc = 1-5, see 
Fig- IHt- In this case we see again in >5'x(w) a peak at 2f2 and in Sy{Lu) a peak at fi. However, the peaks now have a 
small broadening, and small amplitude satellite peaks have appeared at neighboring frequencies. This is evidence of 
another kind of long-term behavior, namely quasiperiodic dynamics. We can corroborate this with the corresponding 
Poincare section shown in Fig. |5Ji. It consists now on a closed one dimensional curve, which means that trajectories 
wind around on a torus, never intersecting itself and yet never quite closing, typical of a quasiperiodic orbit. We have 
also looked at the time dependent estimates of the Lyapunov exponent, X{t). We find that X{t) < for finite t, but 
its absolute value tends to zero for long times as l/t, consistent with quasiperiodic behavior. 

Let us now study the last case, corresponding to regime C. This is done for Idc = 1-66 and lac = 4.05 in Fig.|3;. 
We see that there are broad peaks in the spectrum in both directions, Sx{to) and Sy{uj), and that there is a marked 
increase in the power spectra for low frequencies. Moreover, in Fig. |5f we show the corresponding Poincare section 
which consists on successive points jumping from one region of phase space to another and forming a complex curve, 
which does not seem to close on itself. It is rather difficult to decide from this plot if it corresponds to a quasiperiodic 
orbit or to a low dimensional attractor of a weakly chaotic orbit. We have obtained also the time dependent estimate 
of the Lyapunov exponent also for this case. We find that X{t) > for all t, but its magnitude is decreasing with 
time as 1/t as far as we have been able to observe. The fact that X{t) is possitive for finite t means that there is a 
dynamical instability that causes a seemingly chaotic behavior at intermediate times and a large noise as seen in the 
low frequency power spectrum. However, for long times it is very likely that the system will settle in a quasiperiodic 
dynamics with X{t) 0. In any case, the regime C is very different from the regime A, as can be seen by comparing 
the power spectra of Fig. and Fig. I^;. 

Another view of the dynamics can be obtained by looking at the behavior of the Lyapunov exponent and the noise 
in the region, where a step is expected, as a function of Iac/0- We proceed as follows: for a given lac, we compute 
the set of values of Lyapunov exponents A and low frequency longitudinal noise Sq = lim^^^o Sx{uj) that correspond 
to currents Idc in the region of voltage where a step is expected. (We look at Idc values for which Q — e <Vx < il + e, 
we consider e = 0.005). We plot the resulting set of values of A as a function of lac/^ in Fig. 10a and the values of Sq 
as a functions of lac/^ in Fig. 10b. The vertical lines in the plot correspond to the zeroes of Jo{Iac/^)- We see clearly 
that near these values there are windows of dynamical instability where A > and where the noise Sq is large. In 
the regions of phase-locking we find a couple of interesting results that are worth mentioning, (i) The most negative 
value of the Lyapunov exponent occurs in the middle of the phase-locked step and its magnitude is proportional to 
the step width AS*!, as given by Eq. (|31|) . (ii) The largest value of the noise Sq occurs at the edge of the phase- locked 
step; its magnitude is also proportional to the step width AS*!, as given by Eq. H31|l . 




(41) 



10 



D. Results for Large JJA 

We will now consider the quality of the simple 2x2 model as representing the dynamics of large N x N JJAs. 
It is known that collective effects at high currents may come into play. At high currents, the Z2 symmetry of the 
ground state can be broken because a driving current can induce domain walls.?^^i2M2, Simulations of IV curves with 
the RSJ model and free boundary conditions, for / = 1/2 and T = have reported a chaotic regime at / > /c related 
to the motion of domain walls. It has been shown that open boundary conditions nucleate domain walls leading 
to a critical current lower than de analytic value Ic = 0-35 < — 1 at T = 0.^^ Moreover, Ciria and Giovanellaii 
have shown microscopically that different dynamical states are possible for the longitudinal Shapiro steps. Besides 
the checkerboard ground state configuration, other stable solutions with domain-walls are possible. Then, depending 
on dc current value and history, domain walls can appear, which are not permitted in the four plaquette model. 
Therefore, in order to evaluate to what extent the four plaquette model is valid in the transverse ac driven case, we 
have calculated numerically IV curves for N x N arrays, for TV = 8, 16, 32, 64, with the full RSJ model used before 
in Ref. .27..,3Q. We use periodic boundary conditions in both directions in the presence of an external dc current, 
Ida plus a perpendicular ac current, JacSin(rit). We solve the dynamical equations with time step At — O.Itj 
(tj = 27rci?Ar/o/$o) and total integration time tint = 2^^At after a transient tint/2. We calculate IV curves as a 
function of lac and fi, increasing dc current, /j^,, from checkerboard ground state at Idc — and then decreasing dc 
current, /j^, from the phase configuration obtained at high current. We use a dc current step Aide = 0.01 to obtain 
Ic and Aide = 0.0001 to calculate the step width. 

One of the relevant results with the four plaquette model is the dependence of the critical current with lac/^ for 
high frequencies, as shown in Fig. and Fig. We have also calculated Ic as a function of lac/^ for high fl in 
large JJA arrays. In Fig. llOb we show the case for a particular high frequency value in a 32 x 32 array. We see 
that it has the same behavior as observed in the four plaquette model: Ic{Iac,^) > ^c(0, 0), ranges of lac/^ around 
the maxima of Ic where there is hysteresis, and a quadratic increase with lac for lac/^ 1- We compare with the 
analytical results expected for /J, Eq. and Eq. l|H7|l . which are represented by dot-dashed hnes. We see that 
Jj obtanained numerically for a large array is in excellent agreement with the analytical result for the 2x2 system. 
This is quite reasonable, since /J corresponds to the limit of stability of the checkerboard ground state, which is 
well represented by the 2x2 model. On the other hand, the I^ shows some small deviation from the 2x2 result, 
Ic{N X N) < /c(2 X 2). Also the range in lac/^ where there is hysteresis is bigger in a large system. The current /J- 
corresponds to the low current limit of stability of the moving (non-zero voltage) state. In large systems, the moving 
state can have domain walls, as was found in Ref. 29 30, and the presence of domain walls can lead to a lower I^. 

In order to analyze more quantitatively in which lac/^ ranges the collective effects could be more relevant, we focus 
on two cases: case a corresponding to values that do not show hysteresis in the critical current in a small system, 
but are close to the edge of the lac/^ range of hysteresis, and case b corresponding to values that show hysteresis 
in the 2x2 system. For each case we calculate the critical current both increasing and decreasing the dc drive, and 
therefore they correspond to a t, a i, ^ T and & J, in Fig. llOb. We show in the inset of Fig. llOb the critical currents 
obtained for all these cases as a function of system size, N. In case a, corresponding to the non-hysteretic region, we 
see that there is no size effect in a t np to = 64. Also we see that a 1= a t for N < 32, while for iV = 64 we 
find that hysteresis has appeared and a J,< a t- In the hysteretic region, case 6, size dependent critical currents are 
obtained for b i, while 6 "f is size independent. Moreover, the amplitude of the hysteresis, b 1 —b i weakly increases 
with system size. 

In Fig.llOb we show the range of phase locking A5i as a function of lac/^ for the 32 x 32 array. We find that, when 
there is phase-locking, the numerically obtained ASi is very accurately described by the analytical result of Eq. 
for the 2x2 model. In the inset of Fig. llOb we show the size dependence of ASi for the case marked as c in the 
plot (it corresponds to the same lac/^ of case a of Fig. ITUk). There is no appreciable size dependence. As observed 
in the simulations of the 2x2 system, we also find here that the phase- locking is lost near the zeros of Jo{Iac/^) 
due to the presence of dynamical instabilities. Also we observe that the presence of hysteresis in the critical current 
is nearly coincident with the absence of phase locking. We find that with increasing system size these regimes of 
dynamical instability are amplified in their extension both in their Idc dependence and in their range of lac/^ around 
the zeros of Jo{Iac/^)- This means that the dynamical instabilities detected in the four plaquette system can lead to 
an increased spatiotemporal chaos in larger systems where collective effects are important. 

V. CONCLUSIONS 

It is important to point out that there are no trivial connections between vortex dynamics in the fully frustrated 
JJA and that of a commensurate vortex lattice moving in a rectangular pinning potential in a bulk superconductor^ 
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This is so because the London model with vortices interacting through pair potentials apply to JJAs only in the limit 
of very low vortex density, such that / = Ha? /^q <g; liSi The fully frustrated case represents, in this last respect, 
an interesting limit for studying, where the complete phase field, rather than just the positions of vortices, should be 
taken into account to describe the dynamics. 

We have found transverse phase locking steps in fully frustrated JJA. This new type of (fractional) giant phase- 
locking steps presents marked differences with the well known longitudinal fractional giant Shapiro steps. Particularly, 
the presence of the transverse ac force increases the critical depinning current with respect to the case without ac 
drive (or with a longitudinal ac drive). We have analyzed both analytically and numerically the behavior of the 
steps as a function of ac amplitude lac and frequency Vl. For lac/^ ^ 1, the depinning critical current and the 
phase-locked step width ASi for V = ^lH/2e, increase quadratically with lac- For lac/^ > 1 we have found windows 
of lac/^ where depinning is hysteretic and phase locking is destroyed due to dynamical instabilities. The emergence 
of a weakly chaotic behavior at zero temperature, in a system with non capacitive junctions, is another particular 
characteristic of transverse phase locking which is absent in longitudinal phase locking in overdamped JJA. Comparing 
with the behavior of large fully frustrated arrays we have found that transverse phase locking can be well described 
by an effective four plaquette model, and that collective effects become more important close to the regions of 
dynamical instability of the four-plaquette model. Our results could be observed experimentally in JJA. In particular, 
the enhancement of the critical depinning current with a transverse ac drive could be an interesting experimental 
consequence of the phenomena reported here. 
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FIG. 1: (a) Schematic JJA showing driving currents, and the repeated two-junction by two-junction superlattice unit cell in 
the ground state, (b) Ground state "checkerboard" vorticity for / = Ha^/^o = 1/2. Black squares represent plaquettes with 
one vortex, white squares represent plaquettes without vortices. 
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FIG. 2: Square four-plaquette model, being a, /3, 7, d, the gauge-invariant phase differences. 
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FIG. 3: IV curve of a 2 x 2 JJA for lac = 2.3 and Q, = 1, showing transverse phase locking and critical currents. Inset shows 
a detail of the hysteresis around the critical current. 
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FIG. 4: Critical current, /c, of a 2 x 2 JJA vs ac-amplitude and frequency, lac/^ for high frequencies, f2 > 1. The hysteresis of 
the critical current is demonstrated by (a) the critical current, /J, switching from zero to non-zero voltage state, and (b) the 
critical current, 1^, switching from non-zero to zero voltage states. See Figure 3. Markers are results of numerical simulations 
and lines are the corresponding predictions: (a) equation 1211 1: (b) dashed curve is the maximum of equations II21II and II37I I. 
while the solid curve is the minimum of the two. The critical current, /i, is predicted to follow the solid curve. 
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FIG. 5: Phase-locking of a 2 x 2 JJA at 14 = f2. Markers are results of numerical simulations and lines are the corresponding 
predictions of equations 13H and 11321 . (a) Phase-locking range in dc current, (b) Offset of the phase-locked step relative to 
the Ohmic curve. See Figure |3 
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FIG. 7: First integer step width, ASi, vs ac-amplitude and frequency, lac/^ for intermediate frequencies, Q > 0.5 (a) and low 
frequencies, f2 < 0.5 (b). Inset: Comparison with longitudinal ac-drive for = 1.38. 
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FIG. 8: IV curves for = 1.5 and their corresponding Lyapunov exponents, A: (a) Part of IV curve with phase locking step 
for lac ~ 3.0. Inset: detail of IV curve near the critical current. No hysteresis in Ic is observed, (b) A for lac ~ 3.0. (c) IV 
curve for lac ~ 4.05, no phase locking step is observed. Inset: detail of IV near the critical current: hysteresis in Ic- (d) A for 
lac = 4.05. 
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FIG. 9: Voltage power spectra and Poincare sections for f2 = 1.5 in different I^c regimes: A, B and C (see Fig. |S| (a) 
Transverse, Sy, and longitudinal, Sx, power spectra for Idc = 1-5 and lac ~ 3.0. A regime, (b) Idc = 1.66 and lac = 3.0. B 
regime, (c) Idc = 1.66 and lac = 4.05. C regime. (d),(e) and (f) are the corresponding Poincare sections. Power spectrum Sx 
is plotted displaced in the y-axis for clarity. 
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FIG. 10: Lyapunov exponents A and low frequency noise So for currents Idc giving voltages near V = ^ plotted as a function 
of lac/^. Vertical dashed lines correspond to the zeros of Jo{Iac/^)- (a) Lyapunov exponents. Symbol A indicates value of 
A at the center of the phase-locked step. Dot-dashed line: curve proportional to ASi{Iac/ii) as given by Eq. (glj. (b)L ow 
frewquency noise. Symbol o indicates value of So at the center of the phase-locked step. Dot-dashed line: curve proportional 
to ASi (Jac/n) as given by Eq. CT . 
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FIG. 11: Critical currents Ic and step witdths ASi obtained from numerical simulation in large arrays (32 x 32 junctions) as a 
function of /ac/f2- (a) Ic obtained increasing Idc {ll, •) and decreasing Idc {Ic, Dot-dashed lines show the analytical results 
of Eq. (I2H and Eq. (I37II . Inset shows the size dependence of /J and Ic for system of size N x N, corresponding to the cases a 
and b indicated in the plot, (b) Width of the first integer phase-locked step ASi, obtained numerically for a 32 x 32 array: •, 
and analytical result of Eq. 1311 : dot-dashed line. Inset shows size dependence of A5i for the case c shown in the plot. 



